function [ blendTraits ] = blendCombinedTraits( codes, startSites, endSites, maxGap)
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here
numTraits = size(codes, 1);
blendIdx = 1;

codeStart = 1;
while (codeStart <= numTraits)
    code = codes(codeStart, 1);
    
    codeEnd = codeStart;
    while (codeEnd <= numTraits && codes(codeEnd, 1) == code)
            codeEnd = codeEnd + 1;
    end
    codeEnd = codeEnd - 1;
    %fprintf('examining code %d from %d to %d\n', code, codeStart, codeEnd);
    [codeStarts, codeStartWins, codeEnds, codeEndWins] = blendCodeTraits(startSites(codeStart:codeEnd, 1), endSites(codeStart:codeEnd, 1), maxGap);
    numCodeTraits = size(codeStarts, 1);
    blendCodes(blendIdx:blendIdx+numCodeTraits-1, 1) = code;
    blendStarts(blendIdx:blendIdx+numCodeTraits-1, 1) = codeStarts(:, 1);
    blendStartWindows(blendIdx:blendIdx+numCodeTraits-1, 1) = codeStartWins(:, 1);
    blendEnds(blendIdx:blendIdx+numCodeTraits-1, 1) = codeEnds(:, 1);
    blendEndWindows(blendIdx:blendIdx+numCodeTraits-1, 1) = codeEndWins(:, 1);
    blendIdx = blendIdx + numCodeTraits;
    
    codeStart = codeEnd + 1;
end


sorted = [blendCodes, blendStarts, blendStartWindows, blendEnds, blendEndWindows];
sorted = sortrows(sorted, [1, 2, 3, 4, 5]);

blendCodes(:, 1) = sorted(:, 1);
blendStarts(:, 1) = sorted(:, 2);
blendStartWindows(:, 1) = sorted(:, 3);
blendEnds(:, 1) = sorted(:, 4);
blendEndWindows(:, 1) = sorted(:, 5);
blendTraits = [blendCodes, blendStarts, blendStartWindows, blendEnds, blendEndWindows];

end

function [blendStarts, blendStartWindows, blendEnds, blendEndWindows] = blendCodeTraits(starts, ends, maxGap)
    %maxGap = 3;
    blendIdx = 1;
    lastStart = starts(1, 1);
    blendedByStartStarts(1, 1) = starts(1, 1);
    blendedByStartEnds(1, 1) = ends(1, 1);
    blendedByStartIdx = 1;
    for i = 2:size(starts, 1)
        chrBreak = 0;
%         for j = 1:size(chrEndSites)
%             if (starts(i, 1) > chrEndSites(j, 1) && lastStart <= chrEndSites(j, 1))
%                 chrBreak = 1;
%             end
%         end
        if (chrBreak == 1 || starts(i, 1) - lastStart > maxGap)
            
            [thisBlendStarts, thisBlendStartWins, thisBlendEnds, thisBlendEndWins] = blendByEnd(blendedByStartStarts, blendedByStartEnds);
            numTraits = size(thisBlendStarts, 1);
            blendStarts(blendIdx:blendIdx+numTraits-1, 1) = thisBlendStarts(:, 1);
            blendEnds(blendIdx:blendIdx+numTraits-1, 1) = thisBlendEnds(:, 1);
            blendStartWindows(blendIdx:blendIdx+numTraits-1, 1) = thisBlendStartWins(:, 1);
            blendEndWindows(blendIdx:blendIdx+numTraits-1, 1) = thisBlendEndWins(:, 1);
            blendIdx = blendIdx + numTraits;
            
            blendedByStartStarts = [];
            blendedByStartStarts(1, 1) = starts(i, 1);
            blendedByStartEnds = [];
            blendedByStartEnds(1, 1) = ends(i, 1);
            blendedByStartIdx = 1;
        else
            blendedByStartIdx = blendedByStartIdx + 1;
            blendedByStartStarts(blendedByStartIdx, 1) = starts(i, 1);
            blendedByStartEnds(blendedByStartIdx, 1) = ends(i, 1);
        end
        lastStart = starts(i, 1);
    end
    [thisBlendStarts, thisBlendStartWins, thisBlendEnds, thisBlendEndWins] = blendByEnd(blendedByStartStarts, blendedByStartEnds);
    numTraits = size(thisBlendStarts, 1);
    blendStarts(blendIdx:blendIdx+numTraits-1, 1) = thisBlendStarts(:, 1);
    blendEnds(blendIdx:blendIdx+numTraits-1, 1) = thisBlendEnds(:, 1);
    blendStartWindows(blendIdx:blendIdx+numTraits-1, 1) = thisBlendStartWins(:, 1);
    blendEndWindows(blendIdx:blendIdx+numTraits-1, 1) = thisBlendEndWins(:, 1);
end

function [blendStarts, blendStartWins, blendEnds, blendEndWins] = blendByEnd(starts, ends)
    maxGap = 3;
    comb = [starts, ends];
    comb = sortrows(comb, [2, 1]);
    starts(:, 1) = comb(:, 1);
    ends(:, 1) = comb(:, 2);
    
    blendIdx = 0;
    lastEnd = ends(1, 1);
    blendedByEndStarts(1, 1) = starts(1, 1);
    blendedByEndEnds(1, 1) = ends(1, 1);
    blendedByEndIdx = 1;
    for i = 2:size(ends, 1)
        chrBreak = 0;
%         for j = 1:size(chrEndSites)
%             if (ends(i, 1) > chrEndSites(j, 1) && lastEnd <= chrEndSites(j, 1))
%                 chrBreak = 1;
%             end
%         end
        if (chrBreak == 1 || ends(i, 1) - lastEnd > maxGap)
            blendIdx = blendIdx + 1;
            avgStart = round(mean(blendedByEndStarts));
            avgEnd = round(mean(blendedByEndEnds));
            startWin = max(max(blendedByEndStarts) - avgStart, avgStart - min(blendedByEndStarts));
            endWin = max(max(blendedByEndEnds) - avgEnd, avgEnd - min(blendedByEndEnds));
            blendStarts(blendIdx, 1) = avgStart;
            blendEnds(blendIdx, 1) = avgEnd;
            blendStartWins(blendIdx, 1) = startWin;
            blendEndWins(blendIdx, 1) = endWin;
            
            blendedByEndStarts = [];
            blendedByEndEnds = [];
            blendedByEndStarts(1, 1) = starts(i, 1);
            blendedByEndEnds(1, 1) = ends(i, 1);
            blendedByEndIdx = 1;
        else
            blendedByEndIdx = blendedByEndIdx + 1;
            blendedByEndStarts(blendedByEndIdx, 1) = starts(i, 1);
            blendedByEndEnds(blendedByEndIdx, 1) = ends(i, 1);
        end
        lastEnd = ends(i, 1);
    end
    blendIdx = blendIdx + 1;
    avgStart = round(mean(blendedByEndStarts));
    avgEnd = round(mean(blendedByEndEnds));
    startWin = max(max(blendedByEndStarts) - avgStart, avgStart - min(blendedByEndStarts));
    endWin = max(max(blendedByEndEnds) - avgEnd, avgEnd - min(blendedByEndEnds));
    blendStarts(blendIdx, 1) = avgStart;
    blendEnds(blendIdx, 1) = avgEnd;
    blendStartWins(blendIdx, 1) = startWin;
    blendEndWins(blendIdx, 1) = endWin;
end

